!--------------------------------
! Programa dedicado a calcular
! as 3 raizes de um polinomio de
! terceiro grau.

! Estão implementados 3 metodos:
!
! 1) Chuta erro, verifica:
! Tipo um Hill-Climbing
!
! 2) de Newton-Raphson
! x_{i+1} = x_i - f(x_i)/df(x_i)

! 3) da Secante:
! x_{i+1} = x_i - f(x_i) + ( x_i - x_{i-1} )/( f(x_i) - f(x_{i-1}) )

!
! Programa compilado e testado com:
! $ gfortran -o exerC exerC.f90
! $ ./exerC
! 
! Forma de utilizacao auto-explicativa
! atravez de comentarios.
!
! _o_o_ oOo _o_o_ 28/09/2010


PROGRAM exerC
  IMPLICIT NONE
  REAL (KIND=8), DIMENSION(500,3) :: raizesBD, raizesNR, raizesMS
  REAL (KIND=8), DIMENSION(3) :: raizes, raizesT, FT, FTT, foo
  INTEGER :: N, i
  
  raizes= [-.45, .3, .6]
  raizesBD(1,:) = raizes
  raizesNR(1,:) = raizes
  raizesMS(1,:) = raizes

  ! Lendo os dados
  ! 1) Numero de iteracoes
  write(*,*) "Entre com o numero de numero de iteracoes:"
  read(*,*) N

  i=2
  open(unit=10,file="tab3_out.dat")
  write(10,*) "iter  |  dir1  |  dir2  |  dir3  |  NR1  |  NR2  |  NR3  |  sec1  |  sec2  |  sec3"
  DO WHILE (i<N+2)
      !!!!!!!!!!!!!!!!
      !!!! Metodo da Busca Direta
      ! Fazendo deslocamento aleatorio
      call random_number(foo)
      foo=(foo-.5)*0.1

      ! Aplicando na funcao F
      raizesT=raizesBD(i-1,:)+foo
      FT=abs(  F( raizesBD(i-1,:) )  )
      FTT=abs(  F(raizesT)  )

      raizesBD(i,:) = raizesBD(i-1,:)

      IF ( FTT(1) < FT(1) ) THEN
          raizesBD(i,1)=raizesT(1)
      END IF
      IF ( FTT(2) < FT(2) ) THEN
          raizesBD(i,2)=raizesT(2)
      END IF
      IF ( FTT(3) < FT(3) ) THEN
          raizesBD(i,3)=raizesT(3)
      END IF
      !!!!!!!!!!!!!!!!

      !!!!!!!!!!!!!!!
      ! Metodo de Newton-Raphson
      raizesNR(i,:)=raizesNR(i-1,:) - F(raizesNR(i-1,:))/Fd(raizesNR(i-1,:))
      !!!!!!!!!!!!!!!


      !!!!!!!!!!!!!!
      ! Metodo da secante
      IF (i == 2) THEN
          raizesMS(i,:)=raizesNR(i,:)
      ELSE
          foo=( raizesMS(i-1,:)-raizesMS(i-2,:) )/(  F( raizesMS(i-1,:) )-F( raizesMS(i-2,:) ) )
          raizesMS(i,:)=raizesMS(i-1,:) - F(raizesMS(i-1,:)) * foo
      END IF
      !!!!!!!!!!!!!!

      write(10,*) i-1, raizesBD(i,:), raizesNR(i,:), raizesMS(i,:)
      i=i+1
  END DO
  
  !!!!!!!!!!!!!!!!!!!!!!!
  ! Funcoes utilizadas neste programa
  CONTAINS
  ! Funcao polinomial utilizada verificar convergencia das raizes
  FUNCTION F(x)
    REAL (KIND=8), DIMENSION(3) :: F
    REAL (KIND=8), INTENT(IN), DIMENSION(3) :: x
    F=18*(x**3) -9*(x**2) -5*x +2
  END FUNCTION F

  ! Derivada de F simetrica de 5 pontos
  REAL FUNCTION Fd(x)
    REAL (KIND=8), DIMENSION(3) :: Fd
    REAL (KIND=8), INTENT(IN),DIMENSION(3) :: x
    REAL (KIND=8) :: dx=0.00005
    Fd=-(F(x+2*dx)-8*F(x+dx)+8*F(x-dx)-F(x-2*dx))/(12*dx)
  END FUNCTION Fd
  !!!!!!!!!!!!!!!!!!!!1!!
end PROGRAM